arXiv: 1504.07418v2 [quant-ph] 25 Jan 2016 


Asymptotic properties of the Dirac quantum cellular automaton 

A Perez 

Departament de Ftsica Tedrica and IFIC, 

Universitat de Valencia-CSIC, 

Dr. Moliner 50, 46100-Burjassot, Valencia 
Spain 

We show that the Dirac quantum cellular automaton [Ann. Phys. 354 (2015) 244] shares many 
properties in common with the discrete-time quantum walk. These similarities can be exploited to 
study the automaton as a unitary process that takes place at regular time steps on a one-dimensional 
lattice, in the spirit of general quantum cellular automata. In this way, it becomes an alternative 
to the quantum walk, with a dispersion relation that can be controlled by a parameter, which plays 
a similar role to the coin angle in the quantum walk. The Dirac Hamiltonian is recovered under a 
suitable limit. We provide two independent analytical approximations to the long term probability 
distribution. It is shown that, starting from localized conditions, the asymptotic value of the entropy 
of entanglement between internal and motional degrees of freedom overcomes the known limit that 
is approached by the quantum walk for the same initial conditions, and are similar to the ones 
achieved by highly localized states of the Dirac equation. 
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I. INTRODUCTION 

The connection between physical processes on a lattice, and the corresponding theories in the continuum, is intrigu¬ 
ing and plagued with difficulties and new features [1-3]. Discretization of quantum field theories that are defined on 
the continuum can be regarded as a powerful calculation tool, a paradigmatic example being QCD on a lattice [4], 
that allows for non perturbative calculations, after a suitable extrapolation is made to the limit of vanishing lattice 
spacing. In the case of fermion fields, one encounters problems like the “fermion doubling”, which can be attacked 
in different ways. This clearly shows that the discretization procedure of quantum field theories is not uniquely de¬ 
fined, with different approaches leading to the same limit in the continuum. In particular, this is true for the Dirac 
equation, which describes the relativistic motion of a spin 1/2 particle, and gives rise to interesting phenomena as the 
Zitterhewegung or the Klein paradox [5]. 

A recent paper [6] introduces a Dirac Quantum Cellular Automaton (DQCA), that describes the relativistic dy¬ 
namics of a spin 1/2 particle on a one-dimensional lattice based on some symmetry principles. Quantum cellular 
automata have been studied by several authors (see, for example [7-14]). The model described in [6] can be regarded 
as a particular case of the two component cellular automaton defined in [8], and works as a set of updating rules on 
discrete space-time coordinates, where the time step and lattice spacing are to be identified with the Planck time Tp 
and Planck length Ip, respectively. In the limit of large wavelengths (as compared to Ip) and small masses m mp, 
with mp the Planck mass, the Hamiltonian representing the DQCA approximates the Dirac Hamiltonian. The model 
also accounts for the above mentioned Zitterhewegung and Klein paradox phenomena [15]. 

Also interesting is the fact that the evolution of the probability distribution [6] resembles the one of a discrete 
time Quantum Walk (QW). The QW is the quantum analogue of the classical random walk. As in the case of 
random walks, QWs can appear either under its discrete-time [16] or continuous-time [17] form. Moreover, it has 
been shown that any quantum algorithm can be recast under the form of a QW on a certain graph: QWs can be used 
for universal quantum computation, this being provable for both the continuous [18] and the discrete version [19[. 
Several experimental setups have been already performed to implement the QW [20-33] (for a comprehensive review, 
see [34]). 

In addition to the probability distribution, one immediately finds that the dispersion relation of the QW can be 
mapped into the one corresponding to the DQCA. Last, but not least, both models reproduce the Dirac equation in 
some limit, a property that has been established by several authors in the case of the QW [35-41]. These similarities 
suggest that the two models may share other properties that are worth studying. This is precisely the motivation of 
this paper. We found some subtleties that will be discussed in detail, once the long term evolution has been derived. 
Therefore, we establish a link between a model motivated from a lattice field theory, on the one side, and a process 
(the QW) that plays an important role in the theory of quantum information. 

This paper is organized as follows. In Sect. H we review the general properties of the DQCA and the QW. The 
dispersion relations of both models are discussed in Sect. HI, and we show that the Dirac Hamiltonian is obtained 
from a suitable limit of the DQCA unitary operator. The similarities and differences of the probability distributions 
for both models are analyzed in Sect. IV. In Sect. V, we derive two different approximations to the long term 
probability distribution of the DQCA: We first obtain a simple formula from the r-th moment of the position operator 
at large time steps, which only describes the gross features of the probability distribution, although we can extract 
the correct analytical behavior of the standard deviation. We next obtain an approximate result with the help of the 
stationary phase method, which turns out to work very well, and correctly describes the details of the oscillations in the 
probability. Sect. VI is devoted to the study of the entanglement between the spatial and internal degrees of freedom, 
as quantified by the entropy of entanglement. We will show that, for a localized initial condition, this magnitude 
saturates the allowed maximum value for a two-dimensional Hilbert space, at variance with the lower limiting value 
which is approached by the QW for the same initial conditions. We discuss the similarity of the obtained result with 
highly localized initial states for the Dirac equation. 


II. GENERAL PROPERTIES OF THE DQCA AND THE QW 

A. QW 

The standard QW corresponds to the discrete (both in time and in space) evolution of a one-dimensional quantum 
system (the walker) in a direction which depends on an additional degree of freedom, the chirality, with two possible 
states: “left” |L) or “right” |i?). The global Hilbert space of the system is the tensor product Hs^Hc. Hs is the Hilbert 
space associated to the motion on the line, and it is spanned by the basis {jx = nd) : n € Z}, where d is the lattice 
spacing, usually taken as the unit length. is the chirality (or coin) Hilbert space, defined as a two-dimensional 
space that may correspond, for example, to a spin 1/2 particle, or to a 2-level energy system. Let us call T- (7+) the 
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operators in Hg that move the walker one site to the left (right), and \L){L\, |-R)(i?| the chirality projector operators 
in He- We consider the unitary transformation 

Uqw = {T- 0 \L){L\ + r+ 0 \R){R\} o {/® C(0)} , (1) 

where C{6) is the coin operator, which acts only on the coin space, and I the identity operator in Hg. Any SU(2) 
matrix can be used but, for our purposes, it is sufficient to parametrize 

/cose -sin 9 \ 

^ sint^ COS& J ^ ' 

with 9 e [0,7r/2] a parameter defining the bias of the coin toss. The effect of the unitary operator Uqw on the state 
of the system in one time step r is \tp{t + r)) = UQw\ip{t))- The state vector can be expressed as 

OO 

l^(^)) = X! ® [an{t) |i?) + bn{t) \L)]. (3) 


From the above we obtain 


|■i/'(n, t)) = (nfi|'i/'(t)) = a„(t) |i?) + 5„(t) \L) (4) 

or, in vector notation \ip{n,t)) = {an{t),bn{t))'^■ At any given time step, the probability distribution of the walker 
can be calculated from 


P{n,t) = \an{t)\^ + \bn{t)\^ ■ (5) 

B. DQCA 

As mentioned in the Introduction, the Dirac Quantum Cellular Automaton is an extension of the Dirac field theory 
to the Planck and ultrarelativistic scales [6, 15]. The model is defined by the repeated action of a unitary operator 
UuA that acts on a spinor field with two internal degrees of freedom on a one-dimensional lattice with spacing 

Ip at time intervals rp. Ip and rp being the Planck length and Planck time, respectively. In other words, x = nip, 
t = krp with n,k G Z. The state \ip{t + rp)) of the system at time t -I- rp is related to the one at time t by 

\-tp{t + Tp)) = UDA\'il’{t)). (6) 

If we represent the two internal degrees of freedom by {|i?) ,|L)}, as in the QW, then using similar steps we can define 

l-ipix, t)) = (xl-ipit)) = tppix, t) |i?) -k iPl{x, t) \L), (7) 

where 'ipR{x,t) and 'il)p{x,t) are the right and left spinor components of \'>p(x,t)), respectively. In vector notation, 
ip{x,t) = {i;R{x,t),ipL{x,t))'^. 

Starting from the hypothesis of unitarity, homogeneity of the interaction topology, invariance under time reversal 
and parity, and minimal dimension for a non-identical evolution, one arrives [6] to a unitary operator that can be 
written as: 


Uda = Vi - {T- O \L){L\ + r+ O \R){R\} - (8) 

with f3 a real number. As already proven in [8], in one spatial dimension the conditions of homogeneity and locality 
give rise to a no-go lemma that prevents the existence of nontrivial scalar quantum cellular automata. In other words, 
every band r-diagonal unitary matrix U which commutes with the 1-step translation matrix Tp is also a translation 
matrix Tp for some fc G Z, times a phase. One way to evade the lemma is by combining two consecutive sites of 
the lattice on a cell, and allow the cells to evolve and communicate (the so-called “partitioning/alternating evolution 
rule”). A more interesting possibility is the combination of the two amplitudes of the cell on a field 'tp{x,t) with two 
or more components. This allows to establish connections with field theories (e.g., the Dirac equation). The simplest 
non-trivial case being associated with a value r = 1, and described by a two component field. In this case, under 
the assumptions of locality, unitarity and parity invariance it is shown in the above reference that one arrives to the 
evolution rule 


ip{t -k 1, a;) = x — 1) wo'tpit, x) -k a; -k 1), 


(9) 
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where W-i, wq, and w+i are 2x2 matrices. Any nontrivial solution to the above evolution rule turns out to be 
unitarily equivalent to the choice 


W-i = cos p 


0 i sin 0 \ 

0 cos 9 j ’ 


ic+i = cos p 


cos 0 0 \ 
isin6* 0 j ’ 


Wq = sin p 


sin0 —icos9 
—i cos 9 sin 9 


( 10 ) 


It can be easily shown that the evolution defined by the DQCA Eq. (8) can be obtained from Eqs. (9,10) making 
the choice 9 = 0, P = sinp. Thus, the DQCA can be regarded as a particular case of the r = 1 maps studied in [8]. 
As discussed in this reference, one can also relate this quantum cellular automaton to the one dimensional version of 
Bialynicki-Birula’s unitary cellular automaton for the Dirac equation [7]. 

Similarly to the QW, we can define the spatial probability distribution as 

P{n,t) = \'tjjR{nlp,t)\^ + \tpLinlp,t)\^ . (11) 

We want to establish a connection between both models. To this purpose, we consider the original DQCA as a 
unitary operation taking place on a lattice with arbitrary spacing d at regular time steps r. In other words, we replace 

Ip —y d, Tp —y T. (1^) 

Notice that we depart from the original motivation of the DQCA as a model to describe the relativistic behavior 
of spin 1/2 particle, and consider the automaton as model that can potentially be realized in the laboratory using 
similar setups as for the QW, and constitutes an alternative to the latter. In what follows, we will investigate the 
analogies and differences between both models. 

We notice that the last term in Eq. (8) only acts on the internal degrees of freedom, and does not include any 
displacement on the lattice. As we show later, this introduces some characteristic features on the evolution of the 
DQCA which are at variance with the QW. 


III. DISPERSION RELATION 

Most properties of the QW are better analyzed by switching to the quasi-momentum space [42]. We introduce the 
basis of states {|p) ,p € [— 7r/i/d, 7r/i/d[} defined by 


n= —oo 


(13) 


The unitary operators that govern both the QW and the DQCA become diagonal in this basis. We represent these 
operators by Uqw{p) and Uda{p), respectively. Furthermore, the internal indices can be expressed in the {|i?), \L)} 
basis. With these notations, we obtain 


and 


/ -ipdjh Q \ 

UQwip)=^^ Q Jd/n)Ci9) = 


^-zpd/h Q _f,-^pd|h Q \ 

g*pd/?igin5/ e^pd/^cos9 )' 


Uda{p) 


yfl - p-^e-^Pd/^ -iP 

-iP ^1 - p-^e^pd/f^ 


(14) 


(15) 


In both cases, the eigenvalues can be written as ??+(p) = e where \{p) satisfies the dispersion 

relation 


cos A(p) = cos0cos(pci//i), (16) 

for the QW, and 

cosA(p) = \/l — P"^ cos{pd/h) (17) 

in the case of the DQCA. Therefore, both dispersion relations take the same form, provided that we identify 


cos 


(18) 
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Many features of the time evolution can be obtained directly from the dispersion relation, such as the proportionality 
constant appearing in the asymptotic behavior of the standard deviation [43], or the design of desired asymptotic 
probability distributions in one [44] or more dimensions [45]. Let us notice, however, that even if the correspondence 
defined Eq. (18) is respected, so that both dispersion relations become equivalent, the operators by Uqw{p) and 
UoAip) posses a different structure, which results in the differences that are discussed in the next sections. 

From the unitary operator Eq. (15) one can extract the corresponding Hamiltonian, similarly to [6[. We first write 
t = It, where I G N, and r is the time step. We then define the Hamiltonian H{p) by 


UhAip) = exp[-^ZTiJ(p)]. 

Following this definition, one finds 

H(r)) = f sin{pd/H) _ ^ 

TsinA(p) \ P — ■\/l — /32 sbi{pd/h) 


Let us now rewrite 


P = 


mdc 

h 


(19) 


( 20 ) 


( 21 ) 


with m a parameter with dimensions of mass. The Dirac Hamiltonian is recovered in the limit pd/h <C 1, mdc/h <C 1. 
In this limit, we have sin A(p) ~ A(p), sin{pd/h) ~ pd/h and ~ 1, so that 


Hip) 


d f pc mc^ 
CT V —pc 


( 22 ) 


By choosing the time step and the lattice spacing such that r = d/c, one obtains the Dirac Hamiltonian, where m 
can be identified with the mass of the particle. The latter condition is a reminder of the original model, where d = Ip 
and T = Tp, obviously related by rp = Ip/c. In our proposal, these two parameters are no longer related to the 
Planck scale, but the Dirac dynamics can be recovered within the above restrictions. Considered as a process that 
may approximately simulate a more general wave dynamics via discretization on a lattice, one would need to set the 
length of the discretization, which would determine the correspondence of the parameters of the physical system under 
simulation with the ones of the model (in our case, the value of /3). Also, the value of the time step r is obtained 
from the requirement of a given elapsed time t, and the available number of time steps in the simulation. 


IV. PROBABILITY DISTRIBUTION 

In spite of sharing a common dispersion relation, both models will differ in several other aspects. For our analysis, 
we will fix some of the parameters appearing in these models. The QW will be studied using 0 = 7r/4 in (2), a choice 
that can be mapped to the standard Hadamard coin. Moreover, we adopt the convention d = 1: In this way, we can 
label the site states as |n) in both cases. To allow a comparison with the QW, as discussed in the previous section, 
we will use the value P = 1/^/2 for the following plots. 

We first study the probability distribution, as defined in Eqs. (5) and (11) . Figure 1 shows both probability 
distributions after t = 200 time steps, starting from the initial localized condition 

Hn,0) = Sn,o. (23) 

One observes clear differences: The QW shows its characteristic peaks and a flat distribution in the middle, whereas 
the DQCA features more complicated structures. As it is well known, the probability distribution of the QW vanishes 
at odd (even) sites of the lattice when t is even (odd). This is not true for the DQCA, as observed for t = 10 on the 
same figure, the reason being that the last term in the unitary transformation (8) gives some probability to stay at 
the same position, in contrast to the QW, where the particle is forced to move left and right at each time step. 

Apart from these differences, the figure indicates that both probabilities spread equally with time, at least for large 
time steps. In fact, this is what Fig. (2) shows, where we plot the standard deviation a{t) as a function of t: After a 
few time steps, we obtain the characteristic ballistic a{t) oc t spreading of the QW in both cases. 
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P(n) P(n) P(n) 



Figure 1. (Color online) Left panel: Probability distribution after t = 200 time steps for the QW, where only even sites are 
plotted. The Initial state is localized at the origin, see Eq. (23). Middle panel: Distribution for the DQCA, with all sites 
showed (red dashed line), compared with the long-term approximation, Eq. (33) (black solid line). The right panel shows the 
differences for a smaller {t = 10) time step, where one clearly sees that the probability of the QW vanishes at odd sites of the 
lattice. 



Figure 2. (Color online) Standard deviation <j{t) as a function of the time step t for the QW (blue solid line), and for the 
DQCA (red dashed line). The initial state is localized at the origin, as in the previous figure. 


V. ASYMPTOTIC PROPERTIES 


From the previous section it becomes apparent that, like in the QW case, we can expect well defined properties for 
the DQCA at large time steps. In fact, there are several methods to analytically derive the long term behavior of the 
probability distribution. In what follows, we will consider a localized initial state, as in the previous section. Such 
narrow states in position space would pose problems in the original formulation of the DQCA, i.e. when the model 
is intended to describe the relativistic dynamics of a spin 1/2 particle, if one considers a wave packet width smaller 
than the Compton wavelength of the particle [46]. However, we have recast the DQCA as a discrete time quantum 
process on an ordinary lattice, similar to the QW. In this case, starting from a localized state (as compared to the 
lattice spacing) can be realized in physical implementations. In fact, this the most commonly studied situation, both 
theoretically and experimentally. Therefore, in order to allow for a direct comparison, we consider the localized state 
as our initial state. Also interesting is the study of an initial Gaussian wave packet. As discussed above, this becomes 
a necessity for the original DQCA motivation. The study of Gaussian wave packets for the DQCA has been done in 
[ 6 ]. 
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A. Weak limit 


We first make use of the method developed in [47] (see also [39]) to obtain the convergence of the r-th moment 

E{x^,t) = {m\^^\m), ( 24 ) 

with X the position operator (a different approach, based on combinatorial methods, was used for the QW in [48]). 
Inserting the resolution of the identity in the basis of {|p)} states (13) one obtains 

E{x\t) = [ dp{')p{p,t) I I (25) 

J-TT dp 

where \Yip,t)) = {p \ tpY)) ^ two-component spinor in quasi-momentum space. Using the unitary operator (15) in 
this basis, we can write 

\Y{p, t)) = UYa{p) \Y{p, 0)) • (26) 

The t-th power of Uda{p) is obtained from the spectral theorem: 

UhA{p) = E ■ (27) 

s=±l 


In the latter equation, \{p) is obtained from the dispersion equation (17), and \(j)sip)) t s = ±1 are the two normalized 
eigenvectors of Uda{p), given by [6[: 


\<l>s{p)) 


1 ( a/ I + sv{p) \ 
V2 V - sv{p) ) ’ 


(28) 


where v{p) = ^ = \/l — sin(p)/sinA(p) is the group velocity that arises from (17). Using Eqs. (26) and (27), one 
arrives to the following relation 

{Y{p,t) I {z^Y I YiP,t)) = YY E I (^s(p) I Y(P,0)} Y +o(E-^). ( 29 ) 

In this equation r/'s(p) indicates the derivative with respect to p, and {t)r = t(t — 1) • • • (t — r -I- 1). After dividing by 
U and taking the limit t —> oo, one obtains 

hm E{x^/tYt) = E r dp C-^Y I (Mp) I Y{p,0)) f ■ (30) 

Let us work out the above expression for the initial localized state Eq. (23), for which \Y{p, 0)) = ^ ^ ^ . In this 

case, one finds | {4>s{p) \ YiPi 0)) Y= l/47r, both for s = 1 and s = —1. On the other hand, we can write = sv{p) 

, so that Eq. (30) becomes 

lim E{x^/tYt) = ^ [ dp [vYp) + {-v{p)Y]- ( 31 ) 

t—j-oo 47r 


We next change the integration variable in the first term of the latter expression by inverting the function y = v{p). 
Similarly, we perform the transformation y = —v{p) on the second term. After some algebra, we arrive to the final 
expression 


\im Eix'^/tYt) = - _ dy p{y)yY 

TT 


(32) 


with the notation p' {y) = 


m 


(i-yOVi-/3"-y" 

[— -^/l — /32, -^/l — /?2] with a probability distribution given by 


Eq. (32) implies that the variable x/t is distributed across the interval 


P{y) = ^p'iy) = — - =■ 

TT 7r(l - -/32 - y2 


(33) 



As mentioned above, the same method has been applied to derive the asymptotic probability distribution of the QW. 
We quote this result for comparison. If one starts from a localized state 


■i/'(n,0) = ^ ^ (5„,o, (34) 

with |ap + |5p = 1, and a coin operator as defined in Eq. (2), the corresponding distribution P{y) can be written as 
[39, 48] 


Piy) = 


I sin0| 


[1 _ (|i,p _ |„p _ 


7r(l — y‘^)yJcos^ 9 — ' ' ' cos^ 9 

valid for |y| < | cos6*|. For the particular case Eq. (23), the above formula simplifies to 

I sin0| 


P{y) = 


7r(l — y‘^)yjcos^ 9 — y'^' 


(35) 


(36) 


which coincides with Eq. (33), provided the identification Eq. (18) is made. 

As shown in Fig. 1, our result Eq. (33) provides a simple approximation to the actual probability distribution of the 
DQCA, although it does not reproduce the oscillations seen on the true evolution. However, it can be used to obtain 
the standard deviation at large time steps, as this magnitude does not depend on the details of the distribution. By 
taking r = 2 in Eq. (32) one obtains 


ait) = Wl-\p\, (37) 

which accounts for the ballistic spreading observed in Sect. IV. Indeed, the previous result was expected from the 
similarity of the dispersion relations of the QW (see [43]) and the DQCA, and confirmed by the equivalence of Eqs. 
(33) and (36). 


B. Stationary phase method 


A better approximation to the long-time asymptotic distribution can be obtained following the stationary phase 
method, as used for the QW in [42]. Let us consider a localized initial condition, such that 


|V'(p,o)) = 


1 


V ^ 

and jap + |5p = 1. Making use of Eqs. (26), (27) and (28), we arrive to |V’(p,t)) = (■i/'fl(p, Q, V’i(P) f))^) where 


tpRip, t) = _ [a cos A(p)t — iav{p) sin X{p)t — ib\J\ — v‘^{p) sin A(p)f] 

v 27r 


^hiPit) = -4= [6 cos A(p)t -I- ibvip) sin A(p)t — zoa/I — v‘^{p) sin A(p)t]. 
V 27r 


The corresponding spinor in position space is obtained from 

dp 


tpR,Lin,t)= -^pP^'iljR^Lip,t), 

J —TT V 


where n € Z. Let us introduce the notation a = n/t, and the functions 

^9iip), 

with gi{p) = 1, g 2 {p) = v{p), and gsip) = yjl — v^(p). Then, the above result can be written as 


. -or 27r 


(38) 


(39) 


(40) 


(41) 


(42) 


■0fl(n, t) = oRe {/i(a, t)} — oRe {hict, Q} — iblm{l 3 ia, t)} 


(43) 
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tpL{n,t) = bRe{Ii{a,t)} + bRe{l 2 {a,t)} — ialm{l 3 {a,t)}. 


(44) 


Our goal is to obtain an approximation to the integrals with gi{p) a smooth function. As seen from the 

definition, Eq. (42), the integrand contains an oscillatory exponential, specially for large values of t. For this kind 
of integrals, we can make use of the stationary phase method [49]. Let us consider the phase ^(p, a) = A(p) + ap 
appearing in these integrals, for a given value of a, as a function of p. The basic idea behind the method is to minimize 
these oscillations by expanding the phase ‘f’(p, a) around some convenient point p{a): 


4>(p, a) ~ $(p(a), a) + {p- p{a)) 


9<f>(p, a) 


dp 


p(o 


+ 2^P-pio:)y 


9^4>(p, a) 


dp^ 


p{a) 


(45) 


In the latter equation, it becomes clear that the second term is responsible for strong oscillations, provided that 
the derivative is large. The idea is to minimize these strong oscillations by choosing p{a) such that 

vanishes. Therefore, one needs to look for the roots Pi{a), i = 1,2 ... of the equation 


d^{p, a) 
dp 


v{p) +0 = 0. 


(46) 


Let us first assume a > 0 . After careful inspection, we obtain the roots Pi{a) = —Ps, and P2 {q:) = n +pg, where Ps = 
arccos ^(1 — — q ;^)/[(1 — a^)(l — For the first solution, one needs to replace A(p) —)■ X{ps) = Xs, X (p) —>■ 
A"(ps) = \/l - l 3 ^ — a^{l — a‘^)/\/ 3 \ > 0 , while for the second solution we have to use A(p) —tt + As, A (p) —>■ 
—A (ps). For a < 0 , one needs to change Ps — > —ps- 

As we show below, to our purposes it will suffice to concentrate on Ii{a, t). After substitution of the above results, 
we obtain the following approximation 


/i(Q;,f) ~ 


^27rtA" (ps) 


[e 


2t0(a)+i7r/4 _|_ it{\a\'7T—(p{oc)-\-7r)—i7z 


^4, 


(47) 


where 4>{a) = Ag — |a|ps, the above expression being valid for |q;| < y/l —/3^. To obtain l 2 {a,t) we make use 
of condition (46). It then follows l 2 {a,t) ~ —ali{a,t). Following a similar argument, we arrive to ~ 

yr — a^Ii{a,t). After taking the real part in Eq. (47), and expanding the two resulting cosinus functions, one can 
easily show that Re{Ii{a,t)'\ vanishes whenever t + n is an odd integer number (remember the definition a = n/t, 
where n G Z). In practice, this means that i?e{/i(a, t)} is zero at odd (even) lattice sites when the time step t is even 
(odd). The opposite result is found for Im{Ii{a,t)}: it becomes zero at odd (even) lattice sites when the time step 
t is odd (even). It can be checked that this properties are indeed obeyed by the exact function. In other words, the 
different terms in Eqs. (43,44) “alternate” their contribution, for a given time step, as a function of n, thus obtaining 
a dynamics in the probability distribution that differs from the QW, as already discussed in Sect. IV. 

One might wonder how the above expressions are modified if we want to use a different system of units such that 
d takes an arbitrary value, as defined in Sect. IB. We will not repeat the procedure, and just give the final answer, 
since the above calculations still hold, if one introduces q = pd as the integration variable in (41). In this way, one 
obtains 'iljji^L{nd,t) as the left hand side in Eqs. (43,44). With this modification, the rest of the above results remain 
unchanged (with p replaced by the new variable q). 

In order to test the accuracy of the above approximations, we will analyze the function t), which we extend to 

arbitrary values of x, so that x = nd with n G Z correspond to the lattice sites. A similar argument would apply to the 
extension of the time step t to arbitrary times, by changing t — t/r in the function Ji(J^,t). Notice that changing 
the value of d to d' corresponds to looking for a new value x' in this function, such that x'/d' = xjd. Therefore, it 
is sufficient to consider d = 1 in what follows. Also, within the stationary phase method the functions for 

i = 2,3 are related to in a simple way, so that it suffices for us to consider the latter function. 

Fig. 3 shows the real and imaginary parts of the function for t = 10 time steps, as obtained from direct 

numerical integration of Eq. (42), compared to the obtained approximation Eq. (47). As can be seen from the plots, 
both curves show only an overall resemblance at small number of time steps. However, we show them in order to 
better appreciate the above mentioned parity properties, i.e. in this case the real part vanishes at odd sites, whereas 
the imaginary part does at even sites. The agreement between both functions improves if ones restricts to physical 
sites of the lattice (i.e., for points x = nd such that n G Z). This can be appreciated from Fig. 4, where the same 
results as in Fig. 3 are represented, restricted to lattice sites. On the same figure, one can see that the approximate 
formula for works better for larger time steps, as expected from the stationary phase method, although it 

deviates from the exact value as |q;| approaches the maximum \/l — A similar degree of agreement can be found 
if one extends /i(|,t) to arbitrary (i.e., non integer) values of t. 



10 


Re[/i(x/t,t)] lm[/i(x/t,t)] 




Figure 3. (color online) Real part (left) and imaginary part (right) of the function for t = IQ time steps, as a function 

of X. We adopted the value /3 = \l-\f2 . The dashed-red curve was obtained from a numerical integration of Eq. (42), while 
the blue solid curve corresponds to the obtained approximation Eq. (47). 


Re[;i(n/t,t)] Re[;i(n/t,t)] 



Figure 4. (Color online) Real part of the function where only lattice sites {x = n) are plotted. Blue dots are obtained 

from numerical integration, whereas red squares correspond to the approximate formula. The same value /3 = l/%/2 was used. 
The left panel corresponds to t = 10, while the right panel is for t = 200. In the latter plot, only even sites are represented for 
a better visualization, since the function vanishes at odd sites. 


We have represented in Fig. 5 the probability distribution for the DQCA, as obtained from the stationary phase 
method, compared with the exact evolution, starting from the localized initial condition Eq. (38) with Q = 

b = ^- The plot shows that this approximation works very well, and accurately describes the oscillatory behavior 

of the probability within the limits laj < A detailed analysis shows that the differences in both curves are 

always lower than ^ 5% whenever ta € Z (i.e., for points with support on the lattice). 

The details of these oscillations are better seen for a simple case, corresponding to the initial condition 

V’(n,0) = ^ ^ (5„,o. (48) 

For this particular case, the probability distribution can be expressed, after some algebra, as 

Pia,t) = ^ {(1 + a)^[l + (-l)‘’^*“]cos^(t())(a) -b7r/4) -b (1 - a^)[l - (-1)*+*“] sin^(t(;i(Q;) +7r/4)} . (49) 

This result shows a clear difference with the corresponding result for the QW (c.f. Eq (8) in [42]), where one obtains a 
common factor 1 -b (—!)*+*“, reflecting the parity properties of the QW: for even (odd) t, the probability distribution 
vanishes at odd (even) sites. Instead, Eq. (49) contains a contribution of both even and odd sites at any time step t, 
as detailed above. 
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Figure 5. (Color online) Comparison of the results from the stationary phase method (solid blue line), with the exact evolution 
(dashed red line), after 200 time steps. The initial state is the localized initial condition with o, ~ ^ ~ 


We now return to Fig. 5 to make the following observation. As already observed in Sect. II, the asymptotic form 
of the DQCA looks very similar to the one obtained for the QW. One may wonder about the mathematical reasons 
behind this resemblance. The weak limit showed at the beginning of this Sect, can be used to give an answer, as it 
manages to describe the overall shape of the distribution. As can be seen from Eq. (31), this shape is governed, at 
late times, by the dispersion relation, which can be made to coincide once the relevant parameters (the coin angle for 
the QW, and parameter (3 for the DQCA) are conveniently mapped to each other using (18). Of course, one has to 
remember that details in both distributions are different, as discussed above. 

Similar conclusions can be reached within the stationary phase method. The analysis that was used to obtain 
an approximate expression for the functions Ii{a,t) starts with the expansion Eq. (45), that only depends on the 
dispersion relation. Thus, these functions are the same for the QW and the DQCA, once the mapping (18) is 
established. Now, which precise combinations of the real and imaginary part of these functions one needs is dictated 
by the model (for the DQCA, this combination is given by Eqs. (43,44), whereas for the QW one would need a 
different combination). Again, this explains that we observe a similar shape in the probability distribution, modulo 
the obtained differences. 


VI. ENTANGLEMENT 

A characteristic property of the QW is that entanglement between the coin and spatial degrees of freedom is 
generated as a consequence of the evolution [50“58]. The amount of entanglement is usually quantified using the von 
Neumann entropy of the reduced density matrix of the coin degrees of freedom, after tracing out the spatial ones. 
More precisely, we define this quantity, as a function of the time step t, by 

S{t) = -Tr {pc(t) log 2 Pcit)} , (50) 

where pdt) = the reduced density matrix for the coin space, Tr represents the trace operation 

in this space, and log 2 is the logarithm in base 2. 

As numerically obtained in [50], and proven later in [53], for a Hadamard walk with localized initial conditions the 
asymptotic entanglement is Sum — 0.8720 for all initial coin states, although higher values can be reached by starting 
from non-localized conditions (see also [59]). An obvious question is whether the DQCA is also limited to this amount 
of entanglement, when the evolution starts from the same state. Fig. (6) plots the entropy of entanglement S{t) as 
a function of the time step for both the QW and the DQCA. We immediately see that, for the QW, one approaches 
the predicted value Sum- Interestingly, the DQCA model overcomes this value, and reaches the allowed maximum 
Smax = 1 for a 2-dimensional system, thus indicating that internal and motion degrees of freedom become maximally 
entangled. Such large values of the entanglement are also reached for the one-dimensional Dirac equation with narrow 
initial conditions, for some configurations of the internal degrees of freedom, including the one used in Eq. (23) [36]. 
As already discussed in Sect. V, when considering highly localized states for the Dirac equation, one has to be careful. 
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S(t) 



Figure 6. (Color online) Entropy of entanglement as a function of the time step t for the QW (blue solid line) and for the 
DQCA (red dashed line). The initial state in both cases is localized at the origin, Eq. (23). 


since arbitrarily peaked states are inconsistent with the one-particle approach [46]. One can, however, consistently 
restrict to positive (or negative) energy eigenvalues. For such states, the evolution of highly localized wave packets 
gives rise to maximal entanglement. 

We can get some insight into the above numerical results by obtaining an analytical expression for the reduced 
density matrix Pc{t) in the long term limit. This calculation is more conveniently done in the quasimomentum 
space. We start from the initial state (38), and make use of Eqs. (39,40), which allows us to obtain \'tjj{p,t)) = 
('0k(p, t),'!/’i(p,0)^- Therefore we have 


Pc{t) 



dp\tp{p,t)) (V'(p,i)| • 


(51) 


After expansion, the matrix defined by |■^/;(p, t)) (?/;(p, t)[contains terms of the form sin A(p)t cos A(p)t, sin^ A(p)t, and 
cos^ X{p)t. For large values of t, such terms become highly oscillatory, while the rest of terms that depend on 
the variable p are smooth functions. Thus, we can replace the oscillatory contributions by their averaged value: 
sin\{p)t cos\{p)t —0, sin^ A(p)t —1/2, cos^ A(p)t —1/2. The integral over the resulting expression becomes 
trivial, and we finally obtain 


n =limo(t)=^-(^ 1 ^ 1 ' - 2/3i?e(a6*) A 

2 2pRe{ah*) p \a\^ - {p - 2) \b\^ ) ' 


We can further represent the initial state on the Bloch sphere 


IV'(p, 0 )) 


1 / cos ^ A 

sin^J 


(52) 


(53) 


where 7 G [0, tt] and (p G [0, tt]. Then the above result can be expressed as 


1 / 1 -I- (/3 — 1) cos 7 P cos ip sin 7 

2 I /3cos 93 sin 7 1 —(/3—l)cos 7 


(54) 


As observed from this expression, for the choice 7 = tt/2, together with ip = 7 r/ 2 , 37 r /2 ones has 



(55) 


independently of the parameter p. For such values, then 


lim S(t) = 1 

t—>-oo 


(56) 
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which agrees with the result showed in Fig. ( 6 ), since the state Eq. (23) used for this calculation can be described by 
the values 7 = 7r/2, </? = 7 r/ 2 . 

The differences observed in the amount of entanglement generated within the DQCA, as compared to the QW, 
may have important consequences. The QW has been suggested as a possible device to generate entanglement in 
quantum information processes [60]. On the other hand, the coin can be regarded as a thermodynamic subsystem 
interacting with the lattice. As such, it becomes an interesting scenario to investigate the approach to thermody¬ 
namical equilibrium in quantum systems [61]. We have shown that the DQCA behaves differently to the QW, with 
a dynamics that allows to reach the maximum allowed entanglement. Therefore, it is possible that the transition 
towards equilibrium will show new features. Among these features is the investigation of a non-Markovian behavior 
previous to the asymptotic regime, as already observed for the QW [62]. All these perspectives clearly deserve further 
investigation. 


VII. CONCLUSIONS 


The connection of field theories on a lattice with simpler models that can be used, in some limit, to simulate 
those theories, has proven to be both a useful computational tool, and an avenue towards the understanding of the 
underlying difficulties of the initial system. In this work, we have investigated the time evolution of the Dirac Quantum 
Cellular Automaton, initially proposed as a discretized version that accounts for the motion of a relativistic spin 1/2 
particle in one dimension [6], and compared its properties with those of the Quantum Walk, an important primitive 
for quantum information. We departed from the original motivation of the DQCA, and redefined it as a discrete time 
quantum process taking place on an ordinary lattice, analogously to the QW, that can in principle be implemented 
using similar physical realizations, and allowing to compare both processes on an equal footing. 

The probability distribution looks similar for both systems, with some differences which arise from the fact that 
the DQCA includes a term in the probability amplitude that forces the walker to stay at the original position, at 
variance with the known properties of the QW. In spite of these differences, both probability distributions propagate 
in a similar manner, as clearly shown by the close resemblance of the standard deviation in both cases. Given the 
analogy in the propagation properties, one would expect similar capabilities in applications to quantum algorithms. 
However, in order to reach the full potential of the DQCA, one probably will need to generalize it to more general 
graphs, as in the case of the discrete time QW, in order to look for speedup in quantum search [63, 64], element 
distinctness [65], or even universal quantum computation [19]. 

We have given two analytic approximations to the probability distribution at large time steps. The first one was 
obtained by calculating the generalized momentum of the position operator. In this way, one obtains a simple result 
that only describes the general shape of the distribution, although it suffices to take account for the observed ballistic 
evolution of the standard deviation. On the other hand, the stationary phase method provides a good approximation, 
and clearly shows the effect of the above mentioned “probability to stay” term in the DQCA. 

The analysis of the entanglement between the internal and spatial degrees of freedom reveals that, for some initial 
choices of the coin state localized at one point of the lattice, the DQCA approaches a maximally entangled state . 
This result clearly overcomes the known limiting values of the QW for similar initial conditions. Maximally entangled 
states also appear for some narrow solutions of the Dirac equation in the continuum so that, in this respect, the 
DQCA looks closer to a Dirac particle than the QW. 

To summarize, the DQCA can be regarded as an alternative to the QW, which shares many properties with it, 
while possessing some new distinctive features. At the same time, given its original motivation, it can easily serve 
as a model to illustrate many properties of the Dirac equation, such as the Zitterbewegung and scattering from a 
potential [15]. Of course, an important point is the possibility of experimentally realizing the DQCA. In this respect, 
one might look for setups similar to the ones used to realize the QW, or perhaps make use of some recent ideas on 
quantum dots to implement quantum cellular automata [66-68]. 
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